Design of Oscillator Networks with Enhanced Synchronization Tolerance against Noise 
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Can synchronization properties of a network of identical oscillators in the presence of noise be im- 
proved through appropriate rewiring of its connections? What are the optimal network architectures 
for a given total number of connections? We address these questions by running the optimization 
process, using the stochastic Markov Chain Monte Carlo method with replica exchange to design 
the networks of phase oscillators with the increased tolerance against noise. As we find, the syn- 
chronization of a network, characterized by the Kuramoto order parameter, can be increased up to 
40 %, as compared to that of the randomly generated networks, when the optimization is applied. 
Large ensembles of optimized networks are obtained and their statistical properties are investigated. 
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I. INTRODUCTION 

Synchronization phenomena are ubiquous in various 
fields of science and play an important role in functioning 
of living systems In the last decade, much interest has 
been attracted to studies of complex networks consisting 
of dynamical elements involved in a set of interactions 
0, Q . Particular attention has been paid to problems of 
synchronization in network-organized oscillator systems 
0, [H • Investigations focused on understanding the rela- 
tionship between the topological structure of a network 
and its collective synchronous behavior [3[. Recently, 
synchronization properties of systems formed by phase 
oscillators on static complex networks, such as small- 
world networks [6[ and scale- free networks 0, Hj], have 
been considered. It has also been shown that the ability 
of a network to give rise to synchronous behavior can be 
greatly enhanced by exploiting the topological structure 
emerging from the growth processes [2,LLQ|. However, full 
understanding of how the network topology affects syn- 
chronization of specific dynamical units is still an open 
problem. 

One possible approach is to use evolutionary learn- 
ing mechanisms in order to construct networks with pre- 
scribed dynamical properties. Several models have been 
explored, where dynamical parameters were modified in 
response to the selection pressure via learning algorithms, 
in such a way that the system evolved towards a spec- 
ified goal [Tl -[l3] . This approach can also be employed 
to design phase oscillator networks with desired synchro- 
nization properties. Using heterogeneous oscillators with 
a dispersion of natural frequencies, we have previously 
shown how these elements can be optimally connected, 
by using a given number of links, so that the best syn- 



* Electronic address: yanagita@isc.osakac.ac.jp 



chronization level is achieved [13| . 

Here, our attention is focused on synchronization en- 
hancement in networks of identical phase oscillators in 
the presence of noise. In such systems, noise acting on 
the oscillators competes with the coupling which favors 
the emergence of coherent dynamics [4], LL4[ • The ques- 
tion is how to connect a set of phase oscillators, so that 
the resulting network exhibits the strongest possible syn- 
chronization despite the presence of noise, under the con- 
straint that the total number of available links and, thus, 
the mean connectivity are fixed. 

To design optimal networks, stochastic Markov Chain 
Monte Carlo (MCMC) method with replica exchange 
[l3| is used by us. Large ensembles of optimal networks 
are constructed and their common statistical properties 
are analyzed. As we observe, the typical structure of 
a synchronization-optimized network is strongly depen- 
dent on its connectivity. Sparse optimal networks, with 
a small number of links, tend to display a star-like struc- 
ture. As the connectivity is increased, synchronization- 
optimized networks show a transition to the architectures 
with interlaced cores. 

The paper is organized as follows. In Sec. [Ill we in- 
troduce a model of identical phase oscillators occupying 
nodes of a directionally coupled network and define the 
synchronization measure for this system. The optimiza- 
tion method is also introduced in this section. Construc- 
tion of optimized networks and their statistical analysis 
are performed in Sec.lIIII The results are finally discussed 
in Sec. HV] 



II. THE MODEL AND THE OPTIMIZATION 
METHOD 



For identical oscillators, it is known that, in absence 
of noise, even very weak coupling can lead to complete 
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Below, we consider the effects 
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of noise acting on a network of coupled identical phase 
oscillators, so that the model equations are 
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where are independent white noises, such that 

(&(*)) = and (€i(t)£j(t')) = S 2 SijS(t - £'). Interac- 
tions between the oscillators are specified by the matrix 
w with the elements Wij = 1, if there is a connection, 
and Wij = otherwise. Generally, the connection matrix 
is asymmetric. Note that since the rotation frequencies 
of all oscillators are the same, we can always go into the 
rotational frame Oi \-t Oi — uo^t and thus eliminate the 
term with ujq. Hence, without any loss of generality one 
can set ujq = in Eqs. (pQ). It is known that, for global 
coupling, this model shows a transition to synchroniza- 
tion as the ratio of the coupling strength to the noise 
intensity is increased (see, e.g., [l6[). 

To quantify synchronization of the oscillators, global 
phase 



r(t) 
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will be employed. To measure the degree of synchroniza- 
tion, we numerically integrate Eq. (pQ) with the initial 
conditions 6i(t = 0) = and calculate the average of 
\r(t)\ over a long time T, 



R(w) 



±£\r(t)\dt. (3) 



Our aim is to determine the network w = {wij} which 
would exhibit the highest degree of synchronization, pro- 
vided that the total number K of links is fixed and the 
noise intensity S is given. The network construction can 
be seen as an optimization problem. The optimization 
task is to maximize the order parameter and, possibly, 
bring it to unity by changing the network w. 

To study statistical ensembles of optimized networks, 
the MCMC method [T34HI, which has previously been 
applied to dynamical systems [l3|, l20l - [27| , will be used by 
us. 

We sample networks from the ensemble with the Gibbs 
distribution P(w) ~ exp(/3i?(w)) by the MCMC method. 
To improve the sampling efficiency, we use the Replica 
Exchange Monte Carlo (REMC) algorithm, and the de- 
tails of the algorithm can be seen in [l3| . 

We mainly consider the canonical ensemble average of 
a network function /(-),i.e., 



(Dp = E 



/(w)exp(^(w)) 

Z{fi) 



(4) 



where Z(f3) = exp(/3i?(w)) is the partition function 
and the parameter f3 plays the role of the inverse tem- 
perature. 



III. NUMERICAL INVESTIGATIONS 

To determine the synchronization degree of a given net- 
work at each iteration step of the optimization procedure, 
equations (pQ) were numerically integrated with the time 
increment At = 0.01. Due to limited computational re- 
sources, only relatively small oscillator ensembles of sizes 
TV = 15 are considered in this study. The noise intensity 
is always S = 0.3. 

Initial phases are 0^(0) = 0. Hence, the order param- 
eter at t = is always equal to unity. To construct an 
initial random network with a given number K of connec- 
tions and, thus, with given connectivity p = K/N(N — 1), 
K off-diagonal elements of the matrix w are randomly 
and independently selected and set equal to unity. 

For time averaging, relatively long intervals T = 10000 
were typically used, since the convergence of the order 
parameter is slow. The results did not significantly de- 
pend on T when sufficiently large lengths T were taken. 

In parallel, evolution of M + 1 replicas with different 
inverse temperatures f3 m = 5(3 x m, m = 0, 1, . . . , M has 
been performed (M = 63 and 6/3 = 5). The statisti- 
cal results did not significantly depend on the particular 
choice of inverse temperatures. 

At every five Monte Carlo steps (mcs), the perfor- 
mances of a randomly chosen pair of replicas were com- 
pared and exchanged, as described above. For display 
and statistical analysis, sampling at each every 50 mcs 
after a transient of 5000 mcs has been undertaken. 



A. Optimization at different temperatures 

Synchronization-optimized networks were obtained by 
running evolutionary optimization. In this process, the 
order parameter was progressively increasing until a sat- 
uration state has been reached. Figure Q] gives examples 
of the optimization processes at different temperatures. 
As clearly seen, when using replicas with the larger in- 
verse temperature /?, larger values of the order parame- 
ter could be reached, although the optimization process 
was then slower. This suggests that, for the considered 
problem, the replicas do not actually get trapped in the 
local minima even at large f3 and that already such low- 
temperature replicas can be efficiently used to sample the 
optimized networks. 

After the transients, statistical averaging of the order 
parameter over the ensemble with the Gibbs distribu- 
tion has been performed, according to Eq. ((4]). In Fig. 
[2ja), the averaged order parameter (R)p is displayed as a 
function of the connectivity p for several different inverse 
temperatures j3. The blue solid circle symbols show the 
averaged order parameter corresponding to the replica 
with f3o = 0, i.e. for an infinitely high temperature. We 
see that the averaged order parameter increases with the 
network connectivity p even if the networks are produced 
by only random rewiring. The red open circles show the 
average order parameters for the ensemble corresponding 
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FIG. 1: (Color online). Examples of evolution of the syn- 
chronization order parameter during the optimization process. 
The blue solid, red broken, yellow dotted, green dotted dash 
curves are for the inverse temperatures f3 — /3o, /3i6, ^32 and 
/3m, respectively. The blue solid line (flo = 0) corresponds 
to the networks generated by only random rewiring. The pa- 
rameters are p — 0.1, A — 1.0, 7 — 0.3, M — 63, 5/3 — 5. 



to the replicas with the lowest inverse temperature /3 m- 
Generally, greater order parameters can be obtained by 
running evolution at higher inverse temperatures /?. At 
each connectivity p, the order parameter is gradually in- 
creased with increasing f3 and is approximately saturated 
at /3m- This means that, even if one further increases /?, 
only slight improvements of the averaged order param- 
eter can be expected. Thus, the networks sampled by 
the replica with the largest inverse temperature /3m are 
already yielding a representative optimal ensemble. 

Figure [2fb) shows the ratio (R)p M /(R)/3 of the or- 
der parameters averaged over network ensemble with the 
highest inverse temperature (3m and with the zero in- 
verse temperature (i.e. the ensemble with purely ran- 
dom rewiring) for different connectivities p. Since there 
is no room for the improvement of the order parameter 
when the number of links is small, the ratio tends to 
unity as the connectivity p is decreased. On the other 
hand, when p = 1, global coupling is realized, for which, 
under the chosen coupling strength, full synchronization 
occurs. As evidenced by this Figure, the difference be- 
tween the synchronization capacities of the optimized and 
random networks is most pronounced at the intermediate 
connectivities, for p around 0.1. The noise intensity de- 
pendence for the synchronization capacities is shown in 
Fig[2fc). When the noise intensity is small, the ratio be- 
comes larger and the maximum is shifted to the smaller 
connectivities p. 



B. Collective dynamics 

To analyze differences in the collective dynamics 
of phases oscillators in random and synchronization- 
optimized networks, we have calculated the winding num- 



ber of each oscillator, fli = ^(^(T) — 0^(0)) for many 
realizations of random (sampled by replica with /3q) and 
synchronization-optimized (sampled by replica with /3m) 
networks, and determined the probability distributions 
of winding numbers for both ensembles. As shown in 
Fig. [3j there is a significant difference between these two 
distributions . The probability peak at Qi = for the 
synchronization-optimized ensemble is higher and more 
narrow than that for the random-rewiring ensemble. This 
means that synchronization-optimized networks tend to 
have more elements oscillating with the common fre- 
quency in the presence of the external noises, as com- 
pared with random rewired networks. Thus, elements 
in the synchronization-optimized network behave more 
coherently than those in a random network. 



C. Architectures of synchronization-optimized 
networks 



Several typical synchronization-optimized networks 
are shown in Fig. EJ Their structures strongly depend 
on the number of available connections (the number of 
links is always conserved during an optimization pro- 
cess). When connectivity p is small [Fig. 2] (a)], designed 
networks usually have star structures. The central ele- 
ment acts on a group of periphery elements which have 
no connections among them. Additionally, a number of 
disconnected elements are present. If a larger number of 
links is available [Fig. |4] (b)], a core, formed by a group 
of interconnected elements, becomes formed. There are 
also periphery elements, which are affected by the core, 
but do not influence its dynamics. As the mean con- 
nectivity of the network is increased, the core grows at 
the expense of the periphery elements. Thus, the net- 
work starts to include [Fig. Hl(c)] a relatively large group 
of highly connected elements, with only a few elements 
which are loosely connected and belong to the periphery. 

For a synchronization-optimized network, we have in- 
tegrated the equations (pQ) for a long time, and calculated 
the correlations rji between the phase of a local oscillator 
Oi and that of the global order variable r(t) defined as 
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r(t) exp(i6i)dt 



These quantities show how strongly the dynamics of an 
oscillator i is synchronized with the global signal r(t). 
The nodes in Fig. |4] are colored according to the rescaled 
values rji, i.e., {min^ rji + T^j/lmax^ — min^}. The 
darker color indicates an oscillator having the stronger 
phase correlation with the global signal. 

Figure 2] suggests that the phases of central oscillators 
are strongly correlated with the phase of the global sig- 
nal. In order to check this more clearly, we have divided 
all oscillators into the groups with equal degrees and sep- 
arately determined average correlations with the global 




FIG. 2: (Color online). Average synchronization order parameters (a) and ratios of the order parameters (b) as functions of 
the network connectivity p. The blue filled circles are for the replica /3o, i.e., the ensemble of randomly rewired networks. 
The red squares, yellow diamonds, green triangles, blue inverted triangles, and red open circles are for the replicas with 
P — /3o, /?4, 08, Al6, P32 and /3 m, respectively, (c) Noise intensity dependence of the ratios of the order parameters for /3 m is 
shown for S — 0.2 (blue filled circles), S = 0.3 (red squares), and S — 0.4 (yellow diamonds). Other parameters are the same 
as in Fig. [1] 
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FIG. 3: (Color online). Distributions of the winding num- 
ber for the ensembles of 500 realizations of random rewiring 
networks (sampled by replica with f3o) and of synchronization- 
optimized networks (sampled by replica with /3m). The blue 
circles are for random networks and the red squares are for the 
synchronization-optimized ones. The parameters are same as 
in Fig. [U 

signal for each group. Thus, quantities rjk have been cal- 
culated, 

1 N 
Vk = N Yl VrfkM , 

2^=1 °k,ki i = i 

where hi denotes the total degree of a node z, i.e., 

N N 

hi k~^~ -\- k , k^~ ^ ^ ^hj i ^ y 4 "> 

with k + and k~ being the ingoing and outgoing degrees, 
respectively. 

In Figure [5j phase correlations rjk, averaged over an en- 
semble of synchronization-optimized networks, are plot- 
ted as a function of the degree k for different network 
connectivities p. We see that, on the average, nodes with 



higher degrees are stronger correlated with the global 
signal. Thus, the oscillators having many connections 
act as organizing centers of the synchronization. Fur- 
thermore, as seen in Fig. [5j phase correlations for the 
nodes with the same degree become larger as the connec- 
tivity is increased. This tendency can be understood if 
we take into account that the synchronization-optimized 
networks usually have shallow tree-like structures for the 
smaller connectivities p. As p increases, the network be- 
comes interlaced and has many loops [Fig. H{c)]. Since 
the feedback in a loop enhances the correlation, the av- 
eraged phase correlation of nodes with the same degree 
becomes larger as p increases. 

Note that in a star structure, the central node does 
not receive any signal from other oscillators; thus, the 
phase of the oscillator in the center is only affected by 
the applied noise. On the other hand, when outgoing 
connections from the center to the periphery elements are 
present, the central oscillator effectively acts as a source 
of common noise applied to the peripheral nodes. Re- 
cently, it has been shown that common noise can induce 
synchronization in an ensemble of identical oscillators 
[28L l29| . This phenomenon may be responsible for the 
development of correlations between the peripheral ele- 
ments and the central oscillator. Similar behavior may 
take place when, instead of a single central node, a core 
of highly connected oscillators is present in a network. 



D. Degree distributions 

To statistically investigate architectures of designed 
networks, ingoing and outgoing degrees of their nodes 
have been considered. By sampling over 200 realiza- 
tions from synchronization-optimized ensemble, we have 
obtained the ingoing and outgoing degree distributions 
at p = 0.10, as shown in Fig. [6] For the ensemble 
of random rewiring networks, both ingoing and outgo- 
ing degrees obey the same Poisson distribution (red bro- 
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FIG. 4: (Color online). Examples of synchronization-optimized networks with different connectivities (a) p — 0.05, (b) p = 0.1, 
and (c) p = 0.15. (a'), (b') and (c') correspond to typical random networks with the same connectivity for comparison. The 
nodes are colored according to the phase correlation r/i , and the darker color indicates an oscillator having stronger correlations 
with the global order. The other parameters are same as in Fig. [T] 




6 
k 



10 12 



FIG. 5: (Color online). Correlation between phases of the 
global order and local oscillators as a function of the degree. 
Averaging over 500 realizations sampled by replica with /3m- 
The blue circles and red squares are for p — 0.05 and p — 0.20, 
respectively. The other parameters are same as in Fig. [T] 



ken lines in the figure represent the in- and out-degree 
distributions of networks sampled by the replica with 
f3o). As clearly seen in Fig. [6j most of nodes in the 



synchronization-optimized networks have only one ingo- 
ing connection and no outgoing connections. This indi- 
cates that many periphery nodes exist, consistent with a 
typical realization of synchronization-optimized network 
shown in Fig. |4]^a). Moreover, the outgoing degrees of 
synchronization-optimized networks are distributed more 
broadly than those of random rewiring networks, i.e., 
a long tail in the outgoing connection distribution has 
emerged. This reflects the development of core nodes. 
Hence, there are two principal types of nodes, i.e., core 
and periphery nodes, in the synchronization-optimized 
networks. The core nodes have many outgoing con- 
nections and a smaller number of ingoing connections, 
whereas the periphery nodes tend to have small numbers 
of ingoing connections. 

In order to further investigate the statistics of net- 
work structures as a function of the network connectivity, 
we have calculated the maximum of ingoing and outgo- 
ing degrees of each synchronization-optimized network, 
^max = max^(^ + ) and fc~ ax = max^(/^~), respectively, 
and averaged them over many realizations. In Fig. [7J the 
ratios of averaged maximum ingoing and outgoing de- 
grees of the synchronization-optimized networks to those 
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FIG. 6: (Color online). The distributions of ingoing (a) and outgoing (b) degrees for random and synchronization-optimized 
networks. Each distribution is averaged over 200 network realizations. The blue bars show the distribution for synchronization- 
optimized networks, whereas the red dashed lines are for the random networks. The parameters are same as in Fig. [T] 
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FIG. 7: (Color online). The dependences of relative maxi- 
mum in- and out-going degrees of synchronization-optimized 
networks on their mean connectivity p. The data for out- and 
in-degrees are shown by blue circles and red squares, respec- 
tively. Averaging over 500 realizations of synchronization- 
optimized and random networks. The parameters are the 
same as in Fig. [T] 



of the random networks, i.e., 7+ = (k+ ax ) p M / (k^) p 
and 7- = (&~ ax )/W(^x)/3o, are shown. As p in- 
creases, the ratio of the averaged maximum outgoing de- 
gree of synchronization-optimized networks to that of the 
random-rewiring networks increases steeply and takes the 
maximum in the vicinity of p c = 0.075, while that of the 
outgoing degree (shown by red square symbols) decreases 
and takes the minimum at approximately the same p c . 
In the vicinity of p c , the nodes with a small number of 
ingoing connections and a large number of outgoing con- 
nections (corresponding to the cores) are found in the 
synchronization-optimized networks . 



E. Eigenvalues of the Laplacian matrix 



The Laplacian matrix L for network w is defined as 



u h3 



N 

3=1 



(5) 



Since the considered networks are directed, the eigenval- 
ues of their Laplacian matrices are complex. We can or- 
der the eigenvalues according to the magnitudes of their 
real parts, i.e. as 

= Re(Ai) > Re(A 2 ) > • • • > Re(AAr). 

The eigenvalues of the Laplacian matrix are known to 
play an important role for the synchronizability of oscilla- 
tor networks @- Therefore, we have computed Re A2 and 
Re An/ Re A2 for many realizations of synchronization- 
optimized networks. In Fig. [8] (a), (ReA2)/3 as function 
of the connectivity is shown for different inverse tem- 
peratures /3. It is clearly seen that (ReA2)/3 decreases 
with /3. Since Re A2 determine the inverse relaxation time 
to the synchronized state in oscillator networks [Ej, this 
indicates that the time needed to achieve the synchro- 
nized state decreases with f3. The ratio (Re Xn / Re \2)p 
averaged over the Gibbs ensemble with f3 is shown in 
Fig. [8Kb). The displayed dependencies reveal that the 
ratio, which specifies the to synchronizability, decreases 
as the optimization level, i.e., j3 is increased. 

Recently, fluctuations in the collective signal and os- 
cillation precision were also linked to the eigenvalues of 
the Laplacian matrix [30|, [3l|. The mean intensity of 
the fluctuations of the collective signal in an ensemble 
of components subject to independent Gaussian noises 
can be estimated by the norm of left eigenvector v, cor- 
responding to the zero eigenvalue, vL=0 and normalized 
as YliLi v i ~ 1 ( see details in [31]). When all indepen- 
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FIG. 8: (Color online), (a) ReA 2 , (b) ReA n /ReA 2 , and (c) 
a averaged over the Gibbs ensemble at different inverse tem- 
peratures /3. The circles indicate average over random net- 
works. The squares are averages over the Gibbs ensemble 
with ft = y#8, the diamonds are with ^32 and triangles are 
averages with, i.e., f3 — /3m- The other parameters are same 
as in Fig. [T] 



dent Gaussian noises have the same strength, the mean- 
square dispersion of the collective signal can be estimated 



as a 



f ^2iLi v i- We have computed this property for 
the ensembles of our designed networks. In Fig. |5fc), we 
have shown a as a function of the network connectiv- 
ity for different optimization levels, i.e., different j3. As 
we see, a decreases with /3, implying that collective fluc- 
tuations are suppressed through the optimization. We 
have also estimated the oscillation precision [30| for our 
synchronized-optimized networks. The results indicate a 
tendency to increase the precision at higher optimization 
levels (the figure does not shown in the figure). 

While our networks have been only optimized with re- 
spect to their synchronization ability in the presence of 



noise, the above analysis clearly shows that the designed 
networks turn out to be also optimized with respect to a 
number of other properties. For the designed networks, 
the time of relaxation to the synchronized state in the ab- 
sence of noise is shorter. In the presence of weak noise, 
such networks have lower intensity of fluctuations in the 
collective signal and higher oscillation precision. 



IV. CONCLUSIONS 



We have designed synchronization-optimized networks 
with a fixed number of links for a population of identical 
oscillators under action of independent external noises. 
This has been done by using the Markov Chain Stochastic 
Monte Carlo method complemented by the Replica Ex- 
change algorithm. Large ensembles of networks with im- 
proved synchronization properties have been constructed 
at different mean connectivities and their statistical prop- 
erties have been analyzed by using various characteriza- 
tion tools. 

Our analysis reveals that the architectures leading to 
the improved synchronization of identical oscillators in 
the presence of noise are essentially different from the op- 
timal synchronization architectures for heterogeneous os- 
cillator populations without noise, which have previously 
been studied [l3j|. When the number of available links is 
small, synchronization-optimized networks are typically 
star-shaped structures. As the number of links grows, the 
designed networks are seen to develop dense cores, which 
replace a single central element in the star networks. The 
core expands as the number of available links is increased, 
and eventually the network becomes strongly interlaced. 
The star and core-periphery structures of the designed 
networks can be qualitatively understood, if one takes 
into account that the central elements in such networks 
are effectively operating as the source of common noise 
for the periphery elements. It is known [28|,[29| that com- 
mon noise can induce synchronization in the populations 
of disconnected oscillators or, in our case, in the group 
of periphery elements all connected to the same central 
elements or a central core. 

Thus, we have shown that efficient design of oscilla- 
tor networks with the improved synchronization proper- 
ties is possible. The architectures of such optimal net- 
works strongly depend on the constraints, such as the 
total number of links available. Through the appropriate 
rewiring of a network, a strong gain in the synchroniza- 
tion signal can be achieved. Although our study has been 
performed for a simple system of phase oscillators, sim- 
ilar evolutionary optimization methods can be applied 
to construct networks of different origins, where the dy- 
namics of individual oscillators may be significantly more 
complex. 
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